% THIS CODE ADAPTS NAKAMURA-STEINSSON (NS) PARTIAL EQM MENU COST MODEL 
% TO DESCRIBE FEATURES OF THE PRICING SCHEDULE OF FIRMS UNDER DIFFERENT UNDERLYING EXPECTED/TREND INFLATION

% Solves a simple partial equilibrium menu cost model by Bellman function iteration on a grid.
%
% 1. Constant aggregate demand
% 2. I.i.d. growth of the price level (non-zero mean)
% 3. AR(1)-normal idiosyncratic productivity shocks
% 4. Homoskedastic errors
%
% This model has two state variables: the idiosyncratic level of
% productivity and the firms own relative price.
%
% Jon Steinsson and Emi Nakamura, April 2006
%**************************************************************

cd '/home/user/h112650/paper/replytocgr/MenuCostModelCode/ForDistribution/Partial Equilibrium';

clear;
tic

% Set Parameters - quarterly model (NS use MONTHLY MODEL, I HAVE QUARTERLY SURVEY DATA)
%**************************************************************

% Frequency of price changes
% NS have monthly model and transitions are based on this frequency of adjustment possibilities
% here I want to simulate data also for lower frequencies 
% (proxy for infinite adj cost inbetween periods or rational inattention 
% or other struct features that command lower revision frequency
% therefore for comparability I have to adjust ProbInfl and Proba computed on monthly freq to current frequency of simulation

vecnperiods=[2, 4, 6, 12];


C = 1;      % Aggregate demand
theta = 4;   % Elasticity of demand
vecbeta = 0.96.^(1./vecnperiods);   % Discount factor


% Menu Cost
%%% KK is original NS menu cost
%%%% I consider 3 values


KK = 0.018799*C;
%%vecK = [0.5, 0.6, 0.7, 0.8, 0.9, 1.0, 1.1, 1.2, 1.3, 1.4, 1.5]*KK;
vecK = [ 0.6,  1.0,  1.4]*KK;


% Parameters of the idiosyncratic productivity process (log A_it)
% log(A_it) = rho*log(A_it-1) + eps_it
% where eps_it ~ N(0,sigma_eps^2), i.i.d. 
% The process is truncated at +/- trunc_eps multiples of the unconditional 
% standard deviation of log(A) 
rho = 0.7;
sigma_eps = 0.0425;
trunc_eps = 3.2;

cumrho =[ (1 + rho + rho^2 + rho^3 + rho^4 + rho^5), (1 + rho + rho^2), (1 + rho), 1];

% Number of points on the grid for the two state variables
agridsize = 60;
rpgridsize = 800;


% Max and Min point for the grid of "a"
amax =  trunc_eps*(sigma_eps^2/(1-rho^2))^(1/2);
amin = -trunc_eps*(sigma_eps^2/(1-rho^2))^(1/2);

% rp denotes the real price
% The rp grid is a grid for the log of the real price
rpmax= 0.600;
rpmin=-0.600;
%%% SS RESULTS FOR ALL NPRD, K>=0.6 K<=1.4 AND INFL POPULATE ONLY REL PRICES [-0.3, +0.3]
%%% KEEPING RPSIZE=0.0015 IT MEANS RP GRID SIZE CAN BE 400 PTS INSTEAD OF 800

rpgridsize=400;
rpmax= 0.300;
rpmin=-0.300;

% Create the real price and idiosyncratic productivity vectors
%**************************************************************

a = amin:((amax-amin)/(agridsize-1)):amax; a = a';
rp = rpmin:((rpmax-rpmin)/(rpgridsize-1)):rpmax; rp = rp';

% Calculate Profit function matrix
%**************************************************************
% Prices vary as one moves down each column
% Productivity varies as one moves across each row

rPi = C*(exp(rp).^(-theta+1)*ones(1,agridsize)...
    -exp(rp).^(-theta)*((theta-1)/theta)*(1./exp(a')));


% Create Probability Distribution for the idiosyncratic shock
%**************************************************************
% Proba is a matrix that gives transition probabilities of going
% from state a to state a', i.e. Proba(i,j) = Prob(a' = j | a = i)
% 
% The procedure of Tauchen (1986) is used to approximate the 
% AR(1) process for log(A_it)

MonthlyProba = prob_dist_a_hom(rho,sigma_eps,a);


tmpa=(1:1:rpgridsize)';
tmpb=(1:1:agridsize)';
% fullstate = [RPGRID, AGRID(1)*ONES(SIZE(RPGRID))\ RPGRID, AGRID(2)*ONES(SIZE(RPGRID))\...]
fullstate=[ kron(ones(agridsize,1),tmpa), kron(tmpb,ones(rpgridsize,1))];


vecypi=[0.01, 0.015, 0.020, 0.025, 0.030, 0.035, 0.04, 0.045, 0.05];


nscenarios=size(vecypi,2)*size(vecK,2)*size(vecnperiods,2);

EDP=fullstate;
StatDistr=fullstate;

scencol=0;

for mmm=1:size(vecnperiods,2)
    beta=vecbeta(1,mmm);
    nperiods=vecnperiods(1,mmm)
    
    vecmu=(1+vecypi).^(1/nperiods) - ones(1,size(vecypi,2));
    
    % adjust transitio matrix of tfp shocks
    Proba = MonthlyProba^(12/nperiods);
    asupp =zeros(agridsize,2);
    for ia=1:agridsize
        asupp(ia,1)=find(Proba(ia,:)>0,1,'first');
        asupp(ia,2)=find(Proba(ia,:)>0,1,'last');
    end


for iii=1:size(vecypi,2)

% Parameters of the process followed by the inflation rate (log(P_t/P_t-1))
% log(P_t) - log(P_t-1) = mu + nu_t
% where nu_t ~ N(0,sigma_eta^2), i.i.d.
% this MU is monthly infl rate irrespective of freq of model
% computes monthly transitions for RP
% in simulatiosn will use infl rate adjusted for infra annual freq of price changes
mu = (1+vecypi(iii))^(1/12) - 1; 

% I compute sigma_eta (stdev of monthly inflation) holding sigma/mu = 1.5
% roughly the value implicit in nakamura steinsson original calib:
% monthly infl=0.0021, annual infl=0.0255, monthly sigma = 0.0032,
% cv=0.0032/0.0021

sigma_eta = 1.5*mu;

inflationrate = (1+mu)^nperiods-1;

% Create Probability Distribution for movements in the real price due to inflation
%**************************************************************
% ProbInfl is a matrix that gives the transition probabilities of 
% the real price due to inflation when the nominal price is fixed
%
% ProbInfl gives the probabilities of going from state rp to state 
% rp', i.e. ProbInfl(i,j) = Prob(pr' = j | pr = i)
% 
% The exact distribution is approximated using a 
% method analogous to the Tauchen (1986) method for AR models

% In NS this is MONTHLY TRANSITIONS OF RP
% I adjust it to the current infra annualfrequency

MonthlyProbInfl = prob_dist_infl(mu,sigma_eta,rp);

ProbInfl = MonthlyProbInfl^(12/nperiods);

% here define range of transition conditional on state with positive prob
rpsupp=zeros(rpgridsize,2);
for irp=1:rpgridsize
    rpsupp(irp,1)=find(ProbInfl(irp,:)>0,1,'first');
    rpsupp(irp,2)=find(ProbInfl(irp,:)>0,1,'last');
end

% Initialize ErV matrix (Expected Value next period)
%**************************************************************

ErV = ones(rpgridsize,agridsize)*max(0,mean(mean(rPi))/(1-beta));

tolerance = 0.001;



for k=1:size(vecK,2)
    K = vecK(k);
    
    scencol = scencol + 1;


% Bellman function iteration
%**************************************************************
% Here we iterate ErV

[rV,F,IF] = ValFunIteration(rPi,ErV,beta,K,theta,rp,Proba,ProbInfl,...
    tolerance,1,1);

E1=zeros(rpgridsize,agridsize);
E2=E1; E3=E1; E4=E1; E5=E1; E6=E1; E7=E1; E8=E1; E9=E1; E10=E1; E11=E1; E12=E1;


%% HERE BEGINS COMPUTATION OF QUI INIZIA IL CALCOLOEXPECED VAL OF CHOICE VARIABLE 4 PERIODS AHEAD GIVEN CURRENT STATE SPACE POINT (IRP0, IA0)
if nperiods==12

for irp11 = 1:rpgridsize
    for ia11 = 1:agridsize
        
        for irp12 = rpsupp(irp11,1):rpsupp(irp11,2)
            for ia12 = asupp(ia11,1):asupp(ia11,2)
                
                E12(irp11,ia11) = E12(irp11,ia11) + ProbInfl(irp11,irp12)*Proba(ia11,ia12)*exp(F(irp12,ia12)-rp(irp12));
                
            end
        end
        
    end
end

for irp10=1:rpgridsize
    for ia10=1:agridsize
        
        for irp11 = rpsupp(irp10,1):rpsupp(irp10,2)
            for ia11 = asupp(ia10,1):asupp(ia10,2)
                
                E11(irp10,ia10) = E11(irp10,ia10) + ProbInfl(irp10,irp11)*Proba(ia10,ia11)*exp(F(irp11,ia11)-rp(irp11))*E12(IF(irp11,ia11),ia11);
                
            end
        end
        
    end
end

for irp9=1:rpgridsize
    for ia9=1:agridsize
        
        for irp10 = rpsupp(irp9,1):rpsupp(irp9,2)
            for ia10 = asupp(ia9,1):asupp(ia9,2)
                
                E10(irp9,ia9) = E10(irp9,ia9) + ProbInfl(irp9,irp10)*Proba(ia9,ia10)*exp(F(irp10,ia10)-rp(irp10))*E11(IF(irp10,ia10),ia10);
                
            end
        end
        
    end
end

for irp8=1:rpgridsize
    for ia8=1:agridsize
        
        for irp9 = rpsupp(irp8,1):rpsupp(irp8,2)
            for ia9 = asupp(ia8,1):asupp(ia8,2)
                
                E9(irp8,ia8) = E9(irp8,ia8) + ProbInfl(irp8,irp9)*Proba(ia8,ia9)*exp(F(irp9,ia9)-rp(irp9))*E10(IF(irp9,ia9),ia9);
                
            end
        end
        
    end
end

for irp7=1:rpgridsize
    for ia7=1:agridsize
        
        for irp8 = rpsupp(irp7,1):rpsupp(irp7,2)
            for ia8 = asupp(ia7,1):asupp(ia7,2)
                
                E8(irp7,ia7) = E8(irp7,ia7) + ProbInfl(irp7,irp8)*Proba(ia7,ia8)*exp(F(irp8,ia8)-rp(irp8))*E9(IF(irp8,ia8),ia8);
                
            end
        end
        
    end
end

for irp6=1:rpgridsize
    for ia6=1:agridsize
        
        for irp7 = rpsupp(irp6,1):rpsupp(irp6,2)
            for ia7 = asupp(ia6,1):asupp(ia6,2)
                
                E7(irp6,ia6) = E7(irp6,ia6) + ProbInfl(irp6,irp7)*Proba(ia6,ia7)*exp(F(irp7,ia7)-rp(irp7))*E8(IF(irp7,ia7),ia7);
                
            end
        end
        
    end
end

for irp5=1:rpgridsize
    for ia5=1:agridsize
        
        for irp6 = rpsupp(irp5,1):rpsupp(irp5,2)
            for ia6 = asupp(ia5,1):asupp(ia5,2)
                
                E6(irp5,ia5) = E6(irp5,ia5) + ProbInfl(irp5,irp6)*Proba(ia5,ia6)*exp(F(irp6,ia6)-rp(irp6))*E7(IF(irp6,ia6),ia6);
                
            end
        end
        
    end
end

for irp4=1:rpgridsize
    for ia4=1:agridsize
        
        for irp5 = rpsupp(irp4,1):rpsupp(irp4,2)
            for ia5 = asupp(ia4,1):asupp(ia4,2)
                
                E5(irp4,ia4) = E5(irp4,ia4) + ProbInfl(irp4,irp5)*Proba(ia4,ia5)*exp(F(irp5,ia5)-rp(irp5))*E6(IF(irp5,ia5),ia5);
                
            end
        end
        
    end
end

for irp3 = 1:rpgridsize
    for ia3 = 1:agridsize
        
        for irp4 = rpsupp(irp3,1):rpsupp(irp3,2)
            for ia4 = asupp(ia3,1):asupp(ia3,2)
                
                E4(irp3,ia3) = E4(irp3,ia3) + ProbInfl(irp3,irp4)*Proba(ia3,ia4)*exp(F(irp4,ia4)-rp(irp4))*E5(IF(irp4,ia4),ia4);
                
            end
        end
        
    end
end

for irp2=1:rpgridsize
    for ia2=1:agridsize
        
        for irp3 = rpsupp(irp2,1):rpsupp(irp2,2)
            for ia3 = asupp(ia2,1):asupp(ia2,2)
                
                E3(irp2,ia2) = E3(irp2,ia2) + ProbInfl(irp2,irp3)*Proba(ia2,ia3)*exp(F(irp3,ia3)-rp(irp3))*E4(IF(irp3,ia3),ia3);
                
            end
        end
        
    end
end

for irp1=1:rpgridsize
    for ia1=1:agridsize
        
        for irp2 = rpsupp(irp1,1):rpsupp(irp1,2)
            for ia2 = asupp(ia1,1):asupp(ia1,2)
                
                E2(irp1,ia1) = E2(irp1,ia1) + ProbInfl(irp1,irp2)*Proba(ia1,ia2)*exp(F(irp2,ia2)-rp(irp2))*E3(IF(irp2,ia2),ia2);
                
            end
        end
        
    end
end

for irp0=1:rpgridsize
    for ia0=1:agridsize
        
        for irp1 = rpsupp(irp0,1):rpsupp(irp0,2)
            for ia1 = asupp(ia0,1):asupp(ia0,2)
                % this matrix collectes nperiods ahead expected DP
                
                E1(irp0,ia0) = E1(irp0,ia0) + ProbInfl(irp0,irp1)*Proba(ia0,ia1)*exp(F(irp1,ia1)-rp(irp1))*E2(IF(irp1,ia1),ia1);
                
            end
        end
        
    end
end
 %end if 12 months
 
end

if nperiods==6
   
for irp5=1:rpgridsize
    for ia5=1:agridsize
        
        for irp6 = rpsupp(irp5,1):rpsupp(irp5,2)
            for ia6 = asupp(ia5,1):asupp(ia5,2)
                
                E6(irp5,ia5) = E6(irp5,ia5) + ProbInfl(irp5,irp6)*Proba(ia5,ia6)*exp(F(irp6,ia6)-rp(irp6));
                
            end
        end
        
    end
end

for irp4=1:rpgridsize
    for ia4=1:agridsize
        
        for irp5 = rpsupp(irp4,1):rpsupp(irp4,2)
            for ia5 = asupp(ia4,1):asupp(ia4,2)
                
                E5(irp4,ia4) = E5(irp4,ia4) + ProbInfl(irp4,irp5)*Proba(ia4,ia5)*exp(F(irp5,ia5)-rp(irp5))*E6(IF(irp5,ia5),ia5);
                
            end
        end
        
    end
end

for irp3 = 1:rpgridsize
    for ia3 = 1:agridsize
        
        for irp4 = rpsupp(irp3,1):rpsupp(irp3,2)
            for ia4 = asupp(ia3,1):asupp(ia3,2)
                
                E4(irp3,ia3) = E4(irp3,ia3) + ProbInfl(irp3,irp4)*Proba(ia3,ia4)*exp(F(irp4,ia4)-rp(irp4))*E5(IF(irp4,ia4),ia4);
                
            end
        end
        
    end
end

for irp2=1:rpgridsize
    for ia2=1:agridsize
        
        for irp3 = rpsupp(irp2,1):rpsupp(irp2,2)
            for ia3 = asupp(ia2,1):asupp(ia2,2)
                
                E3(irp2,ia2) = E3(irp2,ia2) + ProbInfl(irp2,irp3)*Proba(ia2,ia3)*exp(F(irp3,ia3)-rp(irp3))*E4(IF(irp3,ia3),ia3);
                
            end
        end
        
    end
end

for irp1=1:rpgridsize
    for ia1=1:agridsize
        
        for irp2 = rpsupp(irp1,1):rpsupp(irp1,2)
            for ia2 = asupp(ia1,1):asupp(ia1,2)
                
                E2(irp1,ia1) = E2(irp1,ia1) + ProbInfl(irp1,irp2)*Proba(ia1,ia2)*exp(F(irp2,ia2)-rp(irp2))*E3(IF(irp2,ia2),ia2);
                
            end
        end
        
    end
end

for irp0=1:rpgridsize
    for ia0=1:agridsize
        
        for irp1 = rpsupp(irp0,1):rpsupp(irp0,2)
            for ia1 = asupp(ia0,1):asupp(ia0,2)
                % this matrix collectes nperiods ahead expected DP
                
                E1(irp0,ia0) = E1(irp0,ia0) + ProbInfl(irp0,irp1)*Proba(ia0,ia1)*exp(F(irp1,ia1)-rp(irp1))*E2(IF(irp1,ia1),ia1);
                
            end
        end
        
    end
end
 
end

if nperiods==4
for irp3 = 1:rpgridsize
    for ia3 = 1:agridsize
        
        for irp4 = rpsupp(irp3,1):rpsupp(irp3,2)
            for ia4 = asupp(ia3,1):asupp(ia3,2)
                
                E4(irp3,ia3) = E4(irp3,ia3) + ProbInfl(irp3,irp4)*Proba(ia3,ia4)*exp(F(irp4,ia4)-rp(irp4));
                
            end
        end
        
    end
end

for irp2=1:rpgridsize
    for ia2=1:agridsize
        
        for irp3 = rpsupp(irp2,1):rpsupp(irp2,2)
            for ia3 = asupp(ia2,1):asupp(ia2,2)
                
                E3(irp2,ia2) = E3(irp2,ia2) + ProbInfl(irp2,irp3)*Proba(ia2,ia3)*exp(F(irp3,ia3)-rp(irp3))*E4(IF(irp3,ia3),ia3);
                
            end
        end
        
    end
end

for irp1=1:rpgridsize
    for ia1=1:agridsize
        
        for irp2 = rpsupp(irp1,1):rpsupp(irp1,2)
            for ia2 = asupp(ia1,1):asupp(ia1,2)
                
                E2(irp1,ia1) = E2(irp1,ia1) + ProbInfl(irp1,irp2)*Proba(ia1,ia2)*exp(F(irp2,ia2)-rp(irp2))*E3(IF(irp2,ia2),ia2);
                
            end
        end
        
    end
end

for irp0=1:rpgridsize
    for ia0=1:agridsize
        
        for irp1 = rpsupp(irp0,1):rpsupp(irp0,2)
            for ia1 = asupp(ia0,1):asupp(ia0,2)
                % this matrix collectes nperiods ahead expected DP
                
                E1(irp0,ia0) = E1(irp0,ia0) + ProbInfl(irp0,irp1)*Proba(ia0,ia1)*exp(F(irp1,ia1)-rp(irp1))*E2(IF(irp1,ia1),ia1);
                
            end
        end
        
    end
end
end

if nperiods==2
for irp1=1:rpgridsize
    for ia1=1:agridsize
        
        for irp2 = rpsupp(irp1,1):rpsupp(irp1,2)
            for ia2 = asupp(ia1,1):asupp(ia1,2)
                
                E2(irp1,ia1) = E2(irp1,ia1) + ProbInfl(irp1,irp2)*Proba(ia1,ia2)*exp(F(irp2,ia2)-rp(irp2));
                
            end
        end
        
    end
end

for irp0=1:rpgridsize
    for ia0=1:agridsize
        
        for irp1 = rpsupp(irp0,1):rpsupp(irp0,2)
            for ia1 = asupp(ia0,1):asupp(ia0,2)
                % this matrix collectes nperiods ahead expected DP
                
                E1(irp0,ia0) = E1(irp0,ia0) + ProbInfl(irp0,irp1)*Proba(ia0,ia1)*exp(F(irp1,ia1)-rp(irp1))*E2(IF(irp1,ia1),ia1);
                
            end
        end
        
    end
end
 
end

if nperiods==1
    for irp0=1:rpgridsize
    for ia0=1:agridsize
        
        for irp1 = rpsupp(irp0,1):rpsupp(irp0,2)
            for ia1 = asupp(ia0,1):asupp(ia0,2)
                % this matrix collectes nperiods ahead expected DP
                
                E1(irp0,ia0) = E1(irp0,ia0) + ProbInfl(irp0,irp1)*Proba(ia0,ia1)*exp(F(irp1,ia1)-rp(irp1));
                
            end
        end
        
    end
end
 
end

    EDP = [EDP, E1(:)];
    
    
% SIMULATE AND COMPUTE DISTR OF REAL PRICES    
NSim = 250000;
burnin = 100;
right_mu=vecmu(iii);
right_rho=rho^(12/nperiods);
right_sigma_eps= (cumrho(mmm)*sigma_eps^2)^0.5;
right_sigma_eta= ((12/nperiods)*(sigma_eta^2))^0.5;


[aSim,InflSim,rpSim,npSim,ncostSim,InflSim_sum] ...
    = DataSimBasic(NSim,burnin,a,rp,right_rho,right_sigma_eps,right_mu,right_sigma_eta,F,1,1);

% Calculate distribution of real price
%****************************************

[SD SD_post] = StaDistNS(rpSim,aSim,InflSim,rp,a);


    StatDistr = [StatDistr, SD(:)];


end    
end
end

save('EDP_MultiK_CV.txt', 'EDP', '-ascii', '-double', '-tabs');        
save('StatDistr_MultiK_CV.txt', 'StatDistr', '-ascii', '-double', '-tabs'); 

scencol

iii

k

mmm


toc

